Why use Mizer?

Species Collapse and Recovery in a Heavily Fished Multispecies System

In this example we will read in a calibrated model for the North Sea, which consists of 12 interacting species and a background resource. All species are predators and prey since they feed according to size-based rules and they encounter each other. Fishing is applied as a logistic curve by length, matched to the maturity length \(l_{50}\) of each species.

More information on how this model was calibrated with historical data can be found in our “mizerHowTo” package (tutorial("HTM3")) and in this blogpost.

First let’s take a look at what the model has used in the historical period (from ICES stock assessments).

The species’ maximum fishing mortality rates (which we call for convenience effort) for Cod and Saithe have declined but are not as low as what would be in line with single-species “Fmsy” (a reference level of fishing deemed sustainable from single-species stock assessements) and the fishing mortality rate for Sprat and Dab is very high.

Then you can take a look at how the modelled species biomasses change through time by running the following:

plotlyBiomass(sim, start_time = 1950, end_time = 2018)

# size spectra averaged for the final 5 years
plotlySpectra(sim,time_range = c(2015:2018),power=2)

Examining changes in the community relative to an unfished state

To be able to examine any changes in the community relative to an unfished state we can re-set effort = 0 and calculate the steady state.

plotlyBiomass(sim0)

We can compare the current size spectra (averaged over 2015-2019) to the unfished size spectra to assess whether there is any evidence of a size-structured trophic cascade due to fishing.

plot_relative_biomass(sim0, sim)

Here we can see the effect of the reduction in large sized individuals of heavily fished species on the other sizes and species in the model, relative to the unfished steady state. The abundance of some (but not all) of the smaller to medium sizes of prey are a lot higher when their larger predators are removed (note the logarithmic scale). Sprat looks to be much lower.

We can do the same with Biomass to see when, if any, of the species collapse. For simplicity, we use < 0.1 of B/B_unfished as a proxy for a reference point for population collapse. We can add other reference points to this kind of plot, for example a simple rule of thumb for B_msy could be 0.5*B_unfished.

Projecting future recovery scenarios

We may now wish to explore the potential recovery of the larger species and sizes in the system. To do this we set up another scenario, where the model starts with the last time step of the fished scenario.

The effort values in 2018 were

Let’s start a new simulation that begins with the effort from 2018 and projects forward for 50 years. We will apply a linear reduction in effort for a selected species to a target value (here assumed for simplicity to be 0.2 with effort expressed in terms of the species’ fishing mortality rate for fully selected sizes).

To do this we need to work with the effort array (time x gear) to enable changes in effort through time, for this scenario. Here we have a ‘gear’ for each species, since effort in this case were annual single-species fishing mortality estimates.

We can see that when we reduce fishing on Sprat it increases the biomass of this species but also affects the biomass of other species in the community.

Are any species collapsed still?

Further Exploration

# allocate species to gear types
gear_params(params)$gear <- c("super_trawler")

params <- setFishing(params, initial_effort = 1)

sim_light <- project(params, effort = .5, t_max = 100)
sim_heavy <- project(params, effort = 1.5, t_max = 100)

plot_relative_biomass(sim_light,sim_heavy)

comment(params@catchability) = NULL # remove comment, needed to add species
params@species_params$biomass_observed = getBiomass(sim)[nrow(sim@n), ] # set biomass_observed

plotCalibration(sim)

plotGrowthCurves(sim, species_panel = TRUE)

species_params <- data.frame(
  species = "Colossal squid",
  w_inf = 3e4,
  w_mat = 2e3,
  beta = 100,
  sigma = 2,
  k_vb = 0.2,
  a = 0.01,
  b = 3
)

new_params = addSpecies(params, species_params)

## No h provided for some species, so using f0 and k_vb to calculate it.

new_params@species_params$biomass_observed[13] = 1e10 # set observed biomass
new_params = steady(new_params)

## Convergence was achieved in 67.5 years.

sim2 = project(new_params, t_max = 5)
plotlyBiomass(sim2, power = 2)

plotBiomassObservedVsModel(sim2)

new_params <- new_params |> calibrateBiomass() |> matchBiomasses() |> steady() |>
  calibrateBiomass() |> matchBiomasses() |> steady() |>
  calibrateBiomass() |> matchBiomasses() |> steady()

## Convergence was achieved in 73.5 years.

## Warning in setBevertonHolt(params): For the following species `erepro` has been
## increased to the smallest possible value: erepro[Colossal squid] = 1.97e-07

## Convergence was achieved in 30 years.

## Warning in setBevertonHolt(params): For the following species `erepro` has been
## increased to the smallest possible value: erepro[Colossal squid] = 1.47e-07

## Convergence was achieved in 4.5 years.

plotBiomassObservedVsModel(new_params)

plotlySpectra(new_params, power = 2)

plotGrowthCurves(new_params, species_panel = TRUE)